Simplified Langevin approach to the Peyrard-Bishop-Dauxois model of DNA 
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A simple Langevin approach is used to study stationary properties of the Peyrard-Bishop-Dauxois 
model for DNA, allowing known properties to be recovered in an easy way. Results are shown for 
the denaturation transition in homogeneous samples, for which some implications, so far overlooked, 
of an analogy with equilibrium wetting transitions are highlighted. This analogy implies that the 
order-parameter, asymptotically, exhibits a second order transition even if it may be very abrupt 
for non-zero values of the stiffness parameter. Not surprisingly, we also find that for heterogeneous 
DNA, within this model the largest bubbles in the pre-melting stage appear in adenine-thymine 
rich regions, while we suggest the possibility of some sort of not strictly local effects owing to the 
merging of bubbles. 
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The DNA thermal melting transition (also called de- 
naturation, coiling, or un-zipping) occurs when, above a 
certain critical temperature, the double-stranded DNA 
molecule unravels into two separate coils, while for 
smaller temperatures (pre-rnelting stage) only localized 
openings or bubbles exist T|. This phase transition is 
of importance for DNA duplication and transcription, 
and many studies have scrutinized its nature (whether 
first or second order), trying to pin down the relevant 
traits of the rich phenomenology experimentally observed 
(a nonexhaustive list of references is 0, 0, 0, H, H, 01)- 
Moreover, it has been suggested that the dynamics of a 
DNA molecule in its pre-melting stage may play a role 
in its own transcription initiation. Indeed, bubbles are 
determined by sequence specificity and they have been 
reported to occur with high probability in the neighbor- 
hood of the, functionally relevant, transcription start site 
(TSS) and near other regulatory sites, facilitating further 
microbiological activity 0, d, Q . 

This relation between thermal dynamics and biological 
functionality has been claimed to be borne out by exper- 
imental data from real promoter DNA sequences and is 
supported by results from a theoretical model (see below) 
[8|, Isl . Even if this might differ from biological, protein 
mediated processes, studies of thermal properties of the 
DNA by itself are a first step forward in understanding 
more complex situations [l[ (see [ifl] for a different view). 

Let us mention some observations in this context, 
which have been the object of recent analyses. Even 
though one would expect that adenine-thymine (AT-)rich 
regions should be more prone to sustain bubbles than 
guanine-cytosine-(GC-)rich ones (as AT pairs bind the 
two strands more weakly than GC ones 1]), counterin- 
tuitive situations in which this is not the case have been 
reported 0, [ll|- In the same vein, the dependence of 
bubble formation on the specific base-pair sequence was 
reported to be highly nonlocal: Upon mutation of two 
AT base-pairs into two (stronger) GC base-pairs near 
the TSS, rendering a specific promoter sequence com- 



pletely inactive for transcription, the opening profiles of 
the original sequence and its mutant variant differed not 
only in the expected suppression of the large thermal 
opening near the TSS, but also in a sizable increase in 
the probability of formation of a bubble at a distant base 
pair [9|. However, subsequent studies using more effi- 
cient methods for the calculation of bubble statistics in 
the Peyrard-Bishop-Dauxois (PBD) model [H, [H did 
not confirm the above non-local scenario, and pointed to 
more localized effects. See [14| . Il5| for recent develop- 
ments on this interesting problem. 

Many of these and other relevant issues have been in- 
vestigated by employing the PBD model [J| (see below). 
The model phenomenology has been profusely analyzed 
by means of various analytical and numerical techniques: 
transfer integral calculations, Monte Carlo simulations, 
molecular dynamics, and Langevin dynamics, and the re- 
sults have been found to pro perly describe experiments 
on the melting transition [ig|, pre-melting bubbles [15J . 
etc. Let us caution that under certain circumstances, 
torsional effects (absent in the PBD model) should be 
included to properly account for some of the described 
phenomenology [Tl [lOl. ITn| . 

In this Brief Report we reconsider the DNA thermal 
denaturation problem analyzing the PBD model [4] by 
means of a different, simplified Langevin approach. This 
strategy allows us to: (i) reproduce numerically in a rel- 
atively easy way the stationary bubble probability dis- 
tribution and other statistical properties for both homo- 
geneous and heterogeneous sequences; (ii) establish an 
analogy with well-known equilibrium wetting problems, 
deeper than previously thought, permitting us to infer 
results about the order of the denaturation transition. 

In the PBD model the stretching of hydrogen-bonds 
between corresponding base-pairs is represented by a set 
of continuous variables {/i„} (at positions n = l,...,iV 
where N is the chain length). The model is defined by 



the following Hamiltoiiian 
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The first term is the kinetic energy for bases of mass 
771. The second one stands for the interaction between 
opposite bases as described by the Morse potential 
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where D„ is the dissociation energy of the 77th base 
pair and a„ denotes the spatial range of the potential. 
Standard, empirically found pair-base-dependent param- 
eter values are customarily employed: D„(AT) — 0.05 
eV, Ai(GC) = 0.075 eV, a„(AT) = 4.2 A-^, and 
a„(GC) = 6.9 A^^ 16]. Finally, the third stacking term 
arises from the interaction between adjacent bases along 
the DNA molecule 3. It reads 

(3) 
where the values of fc, p, and a are determined from 
fittings of experimental DNA denaturation curves [16| : 
k = 0.025 eVA^, p = 2, a ^ 0.35 A'^. The non- 
vanishing stijfness parameter p captures the fact that 
the double-stranded backbone is more rigid than the un- 
wound strands (controlled by a standard elastic inter- 
action). Note that this model includes only transverse 
degrees-of-freedom for nucleotides. 

The average stretching at each site {hn) and its space- 
averaged counterpart {h), as well as {e~^), which can be 
interpreted as the density of closed base-pairs, are the 
standard order-parameters. 

Different scenarios have been reported for the denat- 
uration transition depending on the stiffness parameter 
p and the randomness of the DNA sample. In the sim- 
plest case p — Q ^, the stacking term is harmonic and a 
smooth (second-order) denaturation transition is known 
to occur for both homogeneous and heterogeneous DNA 
[y, ll8[ . On the contrary, non- vanishing p and heteroge- 
neous sequences lead to very abrupt thermal denatura- 
tion curves that exhibit a multistep behavior in line with 
experimental observations [g|. 

The case of nonzero p and homogeneous DNA is still 
unsettled as the transition has been reported to be (i) 
first-order-like yet with a diverging correlation length in 
p^ . Il9| and (ii) second order although very sharp in ap- 
pearance fq]. We shall return to this issue below. Let 
us also remark that, as pointed out in [18|, a continuous 
transition for the order parameter {h) with associated 
critical exponents and a diverging length scale could be 
compatible (if p y^ 0) with the number of bound pairs (77,) 
exhibiting a discontinuity at the transition. 

In evaluating the partition function associated with the 
Hamiltonian Eq.(IT|), the kinetic terms factorize and, as a 
result, can be dropped out if the focus is only on equi- 
librium configurational properties. In such a case, the 



equilibrium state can be recovered from the configura- 
tional part H' of H (including only V and W terms) 
and, therefore, can be reproduced from the stationary 
solution of the associated Langevin equation. 
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where 77 is a Gaussian white noise and a its amplitude. 
In the following, Eq.(j4]) is taken as the starting point for 
study, and an Euler algorithm is used to solve it. This 
differs from previous Langevin studies in that inertial 
terms do not appear, enabling slightly faster computa- 
tional studies. A similar approach was used in [2^. Let 
us stress that the dynamics imposed by Eq.Q is a fic- 
titious one, not related to real DNA dynamics (which is 
not purely relaxational), but leads to the same stationary 
probability distribution as the original one. 

Homogeneous DNA. We begin by studying the case of 
homogeneous samples with only GC base pairs. The tem- 
perature T is the control parameter, and the value of a 
is obtained from the fluctuation-dissipation relation. We 
have run simulations in systems of size 2^^, initializing all 
the base-pairs to h{t = 0) = 2 and letting them evolve 
until a stationary state is reached, (h) was monitored as 
a function of time for zero and nonzero values of p. At low 
temperatures (h) saturates to a finite value whereas at 
high enough temperatures it diverges as t^^^ (see below), 
signaling a phase transition. While for p = a smooth 
(continuous) transition is observed, for p = 2 it is rather 
abrupt (results not shown), being apparently first order. 
The same picture, in line with previous numerical results 
[J], can also be drawn by monitoring {e"'^), but our re- 
sults are not fully conclusive. 

As originally argued in |6|, as hn ~ hn^i, the expo- 
nential factor in Eq.Q can be approximated by e""'*" 
without provoking any significant effect. If p = 0, H' 
is readily recognized (apart from constant terms) as a 
discretized version of the continuous Hamiltonian 
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where 7«i, W2, and k are generic parameters. H^w is the 
standard interfacial Hamiltonian for equilibrium critical 
wetting transitions in the presence of short-ranged forces, 
i.e. the unbinding of the interface separating two coexist- 
ing phases from a wall, which occurs upon increasing the 
temperature [2l[. At this point, we recall that in wetting 
phenomena continuum models are valid approximations 
to lattice models as long as T is above the roughening 
temperature T^, which is Tji = Q in d — 1 [d — 2 bulk). 
Although the connection between wetting and DNA 
denaturation has already been recognized (see, for in- 
stance, [a, [20J) some of its consequences have not been 
fully appreciated. For instance, the set of recently re- 
ported [1^ critical exponents characterizing the DNA 
denaturation transition in the homogeneous cas, {h) ^ 
|(5|-^and^- \6\-'' [where (5 = (T-rc)/Tj, with /3 = -1, 



^ the correlation length, and v = 2, are nothing but the 
two-dimensional critical wetting exponents dating back 
to the early 1980s [2l| . Furthermore, the density of closed 
base pairs scales as {h~^) ~ |<5| (see 0), as corresponds 
to the surface order parameter in a wetting context [2l| . 
Additionally, since in equilibrium wetting the dynamic 
critical exponent z, defined by ^ ~ i^/^, is z = 2, the 
thickness of the wetting layer grows as t^/'^ [24I, in agree- 
ment with the value reported above for the PBD model. 
To the best of our knowledge, these correspondences have 
not been established before. 

More interestingly, the implications of the wetting 
analogy can be extended to the nonzero-p case. In the 
wetting context, a long-standing problem, regarding the 
order of the transition in three-dimensional systems, has 
been recently solved (231] . The original renormalization- 
group calculations led to the prediction of non-universal 
results in blatant disagreement with computational stud- 
ies [2J] and experiments ^] , both of which yield a mean- 
field-like second-order phase transition. An early at- 
tempt to reconcile theory and experiments questioned 
the validity of the effective Hamiltonian Eq.([5]) to de- 
scribe equilibrium wetting and concluded that k in Eq. ([5]) 
should be replaced by a position-dependent sti ffne ss co- 
efficient k{h) = k + w[e-"'' + w'^ahe-^"'' + ■■■ [2^. Cu- 
riously enough, with only the leading correction included 
in fc(/i), this Hamiltonian is the continuous counterpart 
of the PBD one. 

In critical wetting the parameter w'l vanishes at 
the transition point and, according to a linear 
renormalization-group study, only the term proportional 
to W2 is capable of destabilizing the critical wetting tran- 
sition, driving the transition weakly first-order in d = 3 
[2a |. A subsequent investigation allowed the analysis to 
be extended, with the conclusion that a first-order tran- 
sition can appear only for dimensions d ^ 2.41 [27[. Re- 
markably, it has been shown [23| that by including the 
whole series expansion the experimental and computa- 
tional results can be finally reproduced. 

These results can be adapted for homogeneous DNA 
melting. Indeed, by switching on a nonvanishing w'l and 
truncating the series to first order, we do not expect 
the above conclusions to change qualitatively, since it is 
naively expected that w'l plays a similar role to w'2 (the 
detailed proof of this is not straightforward and will be 
published elsewhere). Therefore, using the wetting anal- 
ogy, the one-dimensional melting transition for homoge- 
neous DNA sequences should be asymptotically contin- 
uous for (h), in agreement with some previous transfer 
integral analyses fe] , but in partial disagreement with 
other calculations [ISl, [l^ . Reconciling all these results 
remains an open challenging task. 

Our conclusion about the order of the transition might 
change if we consider versions of the PBD model embed- 
ded in a three-dimensional space [l7[ where bubble en- 
tropic effects are expected to play a crucial role Q . Note 
also that for such three-dimensional models the analogy 
with wetting problems breaks down. 



Heterogeneous DNA. Following the recent literature, 
we have simulated our model for two particular sequences 
of 69 base-pairs: the adeno-associated viral P5 (AAVP5) 
romoter and a mutation of it inactive for transcription 
In the mutant sequence two AT bases located near 
the TSS at positions 48 and 49 are replaced by (more 
tightly bound) GC base pairs. In our analyses a bubble 
is defined as a group of adjacent sites that satisfies the 
condition h > 1.5. To avoid finite-size effects, we use 
periodic boundary conditions on lattices of sizes L = 690 
and 6900 consisting of 10 and 100 replicas, respectively, of 
the same AAVP5 sequence. After sufficient ensemble av- 
eraging, indistinguishable long-time results are obtained 
for both sizes. The bubble distributions for the AAVP5 
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FIG. 1; (Color online). Probability of bubble opening as a 
function of position and bubble size for the AAVP5 promoter 
(top panel) and the mutant P5 promoter (bottom panel) at 
T = 310K. Probabilities in each row are normalized to their 
maximum value as in [I2|. The results are very similar to 
those in 1 131. 



sequence and its mutant are shown in Fig[TJ It can be 
seen that the large bubbles forming around the TSS (top 
panel) are suppressed in the mutant sequence (bottom 
panel) in agreement with experimental observations [8|. 
The effect of the mutation is quite local, in line with that 
obtained in [IjI and in contrast to the first claims [g]. 
Observe, also, that bubbles in the DNA sequence form 
more frequently where AT bases are more abundant, as 
naively expected |14l . Ila ] . Situations in which this is not 
the case (like those reported in [ll|) are likely to be phys- 
ically ascribable to torsional effects [3, [l3| . Our conclu- 
sion is that the local bubble-opening probability within 



the PBD model is controlled by the relative density of 
AT base-pairs, in accordance with |14l.[l5|. 

To explore the possibility of having some sort of nonlo- 
cal effect in bubble formation within the present model, 
consider an artificial chain with a GC-rich region separat- 
ing two AT-rich zones (see Fig 12]). Small bubbles formed 
in the two AT-rich regions might eventually merge to- 
gether, bridging across the GC region as illustrated in 
FiglH This can induce the largest possible bubble to be 
centered around a GC-rich zone, and nonstrictly local 
effects could be generated upon introducing mutations. 
Further research is needed to quantify this mechanism 
and to assess if it is capable of inducing nonlocal effects 
by repetition of the above scenario, which has already 
been discussed in the literature in various forms 1281. 
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FIG. 2: (Color online). Bubble merging over a GC region 
from the openings above two small AT regions. 



In summary, the simple Langevin equation ([4]) gives 
relatively quick access to the stationary properties of the 
PBD model for DNA denaturation. It reproduces many 
known results for the homogeneous case, e.g., for p = 
a continuous transition is obtained. Moreover, we have 
pointed out that the (recently obtained) critical expo- 
nents are well known for the wetting problem. The anal- 
ogy with equilibrium critical wetting can be extended 
using very recent developments to the p ^ case, where 
also a continuous transition is predicted (even if it might 
be a very abrupt one p,[lal)- We have also employed the 
Langevin approach to study the bubble statistics in het- 
erogeneous real sequences, confirming the tendency for 
creation of thermal openings around AT-rich regions. Ac- 
cording to our observations mutations modify the statis- 
tics of bubbles only in a local way. However, nonstrictly- 
local effects due to the merging of bubbles could induce 
large openings in locally GC-rich regions. 



It is our hope that this simple Langevin approach will 
be useful to elucidate other aspects of this fascinating 
field. 



We acknowledge financial support from the Spanish 
MEyC-FEDER, Project No. FIS2005-00791 and from 
Junta de Andalucia as group FQM-165. 



[1] 
[2] 
[3] 



[4] 



[5] 



[6] 

[7] 



[8] 

[9] 
[10] 



R.M. Wartell and A.S. Benight, Phys. Rep. 126, 67 
(1985), and references therein. [11] 

D. Poland and H.A. Scheraga, J. Chem. Phys. 45, 1456 
(1966). [12] 

Y. Kafri, D. Mukamel, and L. Peliti, Phys. Rev. Lett. 
85, 4988 (2000); Eur. Phys. J. B. 27, 135 (2002); E. Car- 
Ion, E. Orlandini, and A. Stella, Phys. Rev. Lett. 88, [13] 
198101 (2002); M.S. Causo, B. Coluzzi, and P. Grass- [14] 
berger, Phys. Rev. E 62, 3958 (2000); S. Cocco and R. 
Monasson, Phys. Rev. Lett. 83, 5178 (1999). [15] 
T. Dauxois, M. Peyrard, and A.R. Bishop, Phys. Rev. E, 
47, R44 (1993). See also, M. Peyrard and A.R. Bishop, 
Phys. Rev. Lett. 62, 2755 (1989); M. Peyrard, Nonlin- 
earity 17, Rl (2004). [16] 
L.-H. Tang and H. Chate, Phys. Rev. Lett. 86, 830 
(2001); T. Hwa et al, Proc. Natl. Acad. Sci. USA. 100, [17] 
4411 (2003); V. Ivanov, Y. Zeng, and G. Zocchi, Phys. 
Rev. E 70, 051907 (2004); M.Y. Azbel, Phys. Rev. A 20, [18] 
1671 (1979); Phys. Rev. E 68, 050901(R) (2003). 
D. Cule, and T. Hwa, Phys. Rev. Lett. 79, 2375 (1997). [19] 
C. J. Benham, Proc. Natl. Acad. Sci. USA 90, 2999 
(1993); R. M. Fye and C. J. Benham, Phys. Rev. E 59, 
3408 (1999). [20] 
C. H. Choi et al, Nucleic Acids Res. 32, 1584 (2004); C. [21] 
H. Choi et al., Phys. Rev. Lett. 96, 239801 (2006). 
G. Kalosakas et al, Europhys. Lett. 68, 127 (2004). 
C. J. Benham and R.R.P. Singh, Phys. Rev. Lett. 97, 



059801 (2006). 

U. Dornberger, M. Leijon, and H. Fritzsche, J. Biol. 
Chem. 274, 6957 (1999). 

T.S. van Erp, S. Cuesta-Lopez, J.-G. Hagmann, and M. 
Peyrard, Phys. Rev. Lett. 95, 218104 (2005); ibid. 96, 
239802 (2006); ibid. 97, 059802 (2006). 
Z. Rapti et al, Europhys. Lett. 74, 540 (2006). 
T.S. van Erp, S. Cuesta-Lopez, and M. Peyrard, Eur. 
Phys. J. E 20, 421 (2006). 

S. Ares et al, Phys. Rev. Lett. 94, 035504 (2005); Z. 
Rapti et al., Phys. Rev. E 73, 051902 (2006); B.S. 
Alexandrov et al. Phys. Rev. E 74, 050901(R) (2006); 
S. Ares and G. Kalosakas, Nano Lett. 7, 307 (2007). 
A. Campa and A. Giansanti, Phys. Rev. E. 58, 3585 
(1998). 

M. Barbi, S. Cocco, and M. Peyrard, Phys. Lett. A 253, 
358 (1999). 

N. Theodorakopoulos, T. Dauxois, and M. Peyrard, 
Phys. Rev. Lett. 85, 6 (2000). 

S. Buyukdagh and M. Joyeux, Phys. Rev. E 73, 051910 
(2006). M. Joyeux and S. Buyukdagh, Phys. Rev. E 72, 
051902 (2005). 

S. Ares and A. Sanchez, Eur. Phys. J. B 56, 253 (2007). 
S. Dietrich, in Phase Transitions and Critical Phenom- 
ena, vol. 12, edited by C. Domb and J. Lebowitz (Aca- 
demic Press, New York, 1988); M. Schick, in Liquids 
at Interfaces, Les Houches Summer School Proceedings 



48, ed. by J. Charvolin, J.F. Joanny, and J. Zinn- Justin [25] D. Ross, D. Bonn, and J. Meunier, Nature 400, 737 

(North Holland, Amsterdam, 1990). (1999). 

[22] R. Lipowsky, J. Phys. A: Math. Gen. 18, L585 (1985). [26] A.J. Jin and M.E. Fisher, Phys. Rev. B 48 2642 (1993). 

[23] A.O. Parry, J.M. Romero-Enrique, and A. Lazarides, [27] C.J. Boulter, Phys. Rev. Lett. 79, 1897 (1997). 

Phys. Rev. Lett. 93, 086104 (2004). [28] T. Novotny et al, Europhys. Lett. 77, 48001 (2007). 
[24] K. Binder, D.P. Landau, and D.M. KroU, Phys. Rev. 

Lett. 56, 2272 (1986). 



